Renormalization of tensor-network states 
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We have discussed the tensor-network representation of classical statistical or interacting quan- 
tum lattice models, and given a comprehensive introduction to the numerical methods we recently 
proposed for studying the tensor-network states/models in two dimensions. A second renormal- 
ization scheme is introduced to take into account the environment contribution in the calculation 
of the partition function of classical tensor network models or the expectation values of quantum 
tensor network states. It improves significantly the accuracy of the coarse grained tensor renormal- 
ization group method. In the study of the quantum tensor-network states, we point out that the 
renormalization effect of the environment can be efficiently and accurately described by the bond 
vector. This, combined with the imaginary time evolution of the wavefunction, provides an accu- 
rate projection method to determine the tensor-network wavfunction. It reduces significantly the 
truncation error and enables a tensor-network state with a large bond dimension, which is difficult 
to be accessed by other methods, to be accurately determined. 



I. INTRODUCTION 

The simulation of strongly correlated quantum or sta- 
tistical systems in two or higher dimensions remains a 
great challenge in physics. This has stimulated great in- 
terest on the investigation of so called tensor network 
states or models in recent years. A tensor network 
state is a high-dimensional generalization of the one- 
dimensional matrix-product stated, [l(| studied by the 
density matrix renormalization group (DMRG).(ll| It 
captures accurately the key feature of entanglement in in- 
teracting quantum systems and is believed to be a good 
starting point for studying correlated systems. On the 
other hand, in a classical statistical system with local in- 
teractions, the Boltzmann weight can be expressed as a 
tensor product and all thermodynamic quantities can be 
determined by studying this equivalent tensor network 
model. 

Strongly correlated quantum spin or fermions systems 
pose some of the most intriguing problems in condensed 
matter theory. Theoretical investigation on these prob- 
lems generally starts from some simplified quantum lat- 
tice models, such as the Heisenberg model, which de- 
scribes the interactions among local spins, or the Hub- 
bard model, which describes the dynamics of band elec- 
trons subject to local Coulomb interactions. It is difficult 
to solve these models because the dimension of Hilbert 
space grows exponentially with the system size. Further- 
more, as quantum fluctuations arc strong in these sys- 
tems, there are no obvious small parameters that can be 
used in perturbative expansions. This is the main ob- 
stacle to the resolution of many fundamental problems, 
such as high-T c superconductivity. 

The quantum Monte Carlo and the DMRG are the two 
most commonly used numerical methods in the study of 



interacting quantum lattice models. The quantum Monte 
Carlo provides an accurate numerical tool for studying 
quantum spin models without frustrations, interacting 
boson models, or some special interacting fcrmion mod- 
els. However, in most of fermion systems, such as the 
Hubbard model away from half-filling, or quantum spin 
models with frustrations, such as the Heisenberg model 
on the Kagome lattice, the quantum Monte-Carlo meth- 
ods are hampered by the "minus-sign" problem. 

The DMRG was proposed by S. White to study ground 
states properties in 1992. fill It has been then extended 
to study thermodvnamic [l2l4l5| as well as dynamic [Tsl— 
[22| properties of quantum lattice models. It can be also 
used to study a system with nonlocal interactions, either 
in the momentum space[23j or other basis spac e (23 . An 
extensive review of DMRG is provided in Ref. [25[ • The 
DMRG minimizes the error in the truncation of basis 
states. It yields extraordinary precise results in one di- 
mension. In two dimensions, the DMRG gives reasonably 
good results for small quantum systems. (2(| HtJ However, 
its truncation error increases dramatically and the result 
becomes less reliable as the lattice size gets larger. 

The DMRG is an iteration method. It starts by di- 
viding a system, called a supcrblock, into two blocks. 
Like an experimental measurement, it uses one of the 
blocks as a needle to probe the states in the other block. 
The response function measured by the DMRG is the 
reduced density matrix. The DMRG works well in one 
dimension because it captures the main feature of the en- 
tanglement between the two blocks. The entanglement 
spectra is determined by the eigenvalues of the reduced 
density matrix, A; . The corresponding entanglement en- 
tropy is defined by 



(1) 
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As pointed out in Ref. |2jj , the DMRG is to maximize the 
entanglement entropy during the basis truncation and 
can be regarded as a maximal entropy method. 
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For a gapped quantum system, the entanglement en- 
tropy is believed to satisfy an area law, namely the en- 
tanglement entropy scales with the cross section between 
the system and environment blocks [28j. For a matrix 
product state, the entanglement entropy is bounded by 



S < log A 



(2) 



where D is the matrix dimension. Thus to accurately rep- 
resent a quantum state by a matrix product, the minimal 
dimension of the matrix, above which the wavefunction 
converges exponentially, must grow exponentially with 
the surface cross section. In one dimension, the surface 
contains just two points, a small number of states D is 
sufficient to make DMRG extremely accurate. 

In two or higher dimensions, the bond dimension that 
is needed for describing accurately a ground state grows 
exponentially with the system size. It is difficult to use 
a matrix product wavefunction to represent faithfully a 
quantum state in two or higher dimensions. That is why 
the DMRG can yield reliable results only in small lattice 
systems in two dimensions. It is an intrinsic barrier in the 
application of the DMRG in two or higher dimensions. 

At a critical point, there is a logarithmic correction to 
the area law.[29l.l30j This logarithmic correction does not 
affect much on the calculation of short range correlation 
functions. However, as it diverges with the system size, 
it can affect strongly on the calculation of long range 
correlation functions. The minimal matrix dimension will 
grow exponentially with the system size. This is why the 
DMRG works better in a gapped system, for example 
the spin-1 Hcisenbcrg model, than in a critical system, 
for example the spin-1/2 Heisenberg model. 

In 1995, Ostlund and Rommer pointed out that the 
wavefunction generated by the DMRG iteration is a 
matrix-product state. flol| A matrix product state can be 
also taken as a variational wavefunction. It can be accu- 
rately determined by the DMRG in one dimension. 

The DMRG is a one-dimensional algorithm. To ap- 
ply the DMRG in more than one dimension, one needs 
to map a two- or higher-dimensional lattice into a one- 
dimensional one by introducing nonlocal coupling. [2?) 
However, locality is the key of the great performance of 
the DMRG in one dimension. To capture the entangle- 
ment between any two neighboring sites, one needs to 
keep the locality of interactions. 

A tensor-network state is a high-dimensional extension 
of the one-dimensional matrix-product state. It keeps 
the locality of interactions and captures the main feature 
of the area lawQ, since the number of entangled bonds 
between the systems and environment blocks is propor- 
tional to the interface area. A tensor product representa- 
tion reflects accurately the short-range entanglement of 
a quantum state. Like in a matrix product state in one 
dimension, the minimal bond dimension that is needed 
for accurately describing a quantum state does not de- 
pend on the system size in a gapped system. This is an 
advantage of using the tensor representation. 



The tensor-network representation of quantum states 
is in fact not new. A typical example of tensor-network 
state is the valence-bond-solid (VBS) state that waspro- 
posed by Afflect et. al more than two decades ago[lJ. It 
is the exact ground state of a class of two-dimensional 
quantum antifcrromagnets. A tensor product wavefunc- 
tion was also constructed by Niggemann et al. to study 
the ground state properties of the Heisenberg model on 
the honeycomb lattice [3] . They pointed out that the cal- 
culation of expectation values of tensor product states is 
equivalent to evaluating a classical partition function. A 
more general ansatz of tensor network states was sug- 
gested by Sierra and Martin-Delgado[31|. The use of 
the tensor network state as a variational wavefunction 
for the three-dimensional classical lattice model was sug- 
gested by Nishino.[32[ In 2004, the idea of tensor net- 
work state was discussed by Cirac and coworkers under 
the name of projected entangled-pair states (PEPS)Q. 
Their work has attracted great attention because it re- 
veals more clearly the physical picture embedded in the 
tensor-network representation of quantum states. 

Even before the tensor-network wavefunction was in- 
troduced in quantum systems, a tensor-network repre- 
sentation of the partition function for classical lattice 
models has been widely used in statistical physics. A 
vertex model[33| is such an example. Furthermore, as 
will be discussed in Sec. |ITJ all classical lattice models 
with local interactions can be written as tensor-network 
models. Unlike in a quantum system, the local tensor 
can be readily determined from its Hamiltonian. 

In the past decade, a number of variational approaches 
have been suggested to determine the tensor-product 
wavefunction in quantum systems. [3II However, the 
number of free tensor elements that can be handled by a 
variational approach is too small to resolve any physical 
problems that are difficult to be resolved by other meth- 
ods. This is because in a variational approach, the max- 
imal number of free variables can be accurately and effi- 
ciently determined from the minimization of the expecta- 
tion value of Hamiltonian is generally less than I0 2 ~ 10 3 . 
This has limited the bond dimension to be generally less 
than 5. 

In Ref. Q, we proposed a projection method (or power 
method) to evaluate the tensor-network wavefunction of 
the ground state. In this method, the projection opera- 
tor is applied iteratively to a random initial wavefunction 
by utilizing the Trotter-Suzuki decomposition. The con- 
tribution of the environment is approximately taken into 
account by the diagonal bond matrix. This method is 
quite efficient and accurate, since in the evaluation of 
the wavefunction by the imaginary time evolution, the 
Trotter and truncation errors do not accumulate in the 
iteration. At each step of iteration, a singular value de- 
composition for a matrix of dimension dD 2 on a hon- 
eycomb lattice, or dD 3 on a square or Kogome lattice, 
needs to be done. Here d is the dimension of the local 
basis set and D is the bond dimension of the local tensor. 
For the spin-1/2 Heisenberg model d = 2. Nowadays a 
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matrix of dimension ~ 10 4 can be cfficiciitly diagonal- 
ized by a desktop computer. Therefore, a tensor-network 
state with the bond dimension up to D = 70 on a honey- 
comb lattice, or D = 17 on a square or Kagome lattice, 
can be calculated for the S=l/2 Heisenberg model. 

The above projection method is similar that used by 
Vidal in one dimension [36]. But there is a significant 
difference between one and two dimensions. In one di- 
mension, there is a canonical representation of the matrix 
product state in which the bond vector is preciously the 
singular spectrum of Schmidt decomposition. If the ma- 
trix product state is canonical, then by cutting a bond, 
the left and right matrix product states are strictly or- 
thogonal. The projection method introduced in Ref. [3(j| 
can be naturally imposed by applying the orthogonality 
condition. However, in two dimensions, this kind of or- 
thogonality condition does not exist and the bond vector 
is no longer the singular spectrum of any kind of Schmidt 
decompostion since to separate a point one needs to cut 
at least four bonds on a square lattice. Our approach 
is to take the bond vector as an effective entanglement 
measure of environment. It goes beyond the canonical 
representation of tensor-network wavef unctions. 

To evaluate the expectation values of physical observ- 
ables from a given tensor-network wavefunction, a num- 
ber of approaches, including the transfer matrix renor- 
malization group introduced in Refs. [13I-H5I] or the varia- 
tional Monte Carlo mehtods(35j. can be used. In Ref. 
we adopted the coarse graining tensor rcnormalization 
group (TRG) method proposed by Levin and Nave [J]. 
From the calculation, we find that the TRG can indeed 
produce qualitatively correct results when the bond di- 
mension of tensors is small. However, the truncation er- 
ror in the TRG iteration grows rapidly with the bond 
dimension of local tensors. This leads to a big error in 
the calculation of expectation values. In particular, the 
ground state energy and other physical quantities do not 
converge in the large bond dimension limit. 

In the TRG method of Levin and Nave, the singular- 
value spectra of a local matrix M defined by a product of 
two neighboring local tensors is renormalized in the trun- 
cation of basis space. This provides a local optimization 
of the truncation space since it does not consider the in- 
fluence of the environment tensors which are defined from 
the whole lattice by excluding the two tensors on which 
M is defined. Similar as in the DMRG, the interplay 
between the system (M here) and environment tensors 
can modify dramatically the singular- value spectra of M. 
Thus the renormalization effect of the environment to M 
should be considered to globally optimize the accuracy 
of physical observable. 

In Ref. Q, we introduced a second renormalization 
method of tensor network model (abbreviated as SRG) 
to study this renormalization effect of environment. This 
method can be used to evaluate the partition functions of 
classical statistical models and the expectation values of 
quantum tensor-network wavef unctions. It improves sig- 
nificantly the accuracy of TRG. For example, for the two- 



dimensional Ising model, we found that this method can 
improve the accuracy of the TRG by more than two or- 
ders of magnitude in the vicinity of the critical point and 
more than five orders of magnitude away from the criti- 
cal point by keeping D cut = 24 states. This method also 
provides a powerful numerical tool for accurately evalu- 
ating the expectation values of tensor-network states in 
a quantum lattice model. In contrast to the TRG, we 
find that the SRG results for the ground state energy of 
the Heisenberg model on a honeycomb lattice converge 
quickly with increasing bond dimension D. 

The SRG together with the projection method intro- 
duced in Ref. J5| establishes an efficient numerical tech- 
nique for studying the tensor network states of quantum 
lattice models in two dimensions. The SRG alone also 
provides an accurate numerical tool for studying ther- 
modynamics of classical lattice models. In this paper, 
we will give a comprehensive introduction to these meth- 
ods. 

This paper is arranged as follows. Sec. [H] is devoted 
to a discussion on the tensor-network representation of 
classical lattice models. Two methods are introduced 
to convert a classical lattice model into a tensor-network 
model. One is to define the local tensor in the dual lattice 
by performing a duality transformation, and the other is 
to define the local tensor by singular value decomposition 
in the original lattice. Sec. IIIII starts by a brief review of 
the TRG method. A detailed introduction to the SRG 
is then given. Sec. IIVI describes the iterative projection 
method for determining the tensor-network wavefunction 
of ground state in a quantum system. The calculation 
of the expectation values of tensor-network states is dis- 
cussed in Sec.|V] Sec. [VI] is a summary. 

For reference, we summarize a few notations used in 
this work: d is the dimension of local physical basis set. 
D is the bond dimension of a matrix or tensor product 
state in a quantum system. D c is the bond dimension of 
a classical tensor network model. In a quantum system, 
D c = D 2 is the bond dimension of the tensor used in 
the evaluation of physical expectation values. D cut is the 
bond dimension of the tensor retained in the TRG or 
SRG calculations. 



II. TENSOR NETWORK REPRESENTATION 
OF CLASSICAL STATISTICAL MODELS 

In this section, we discuss about the tensor-network 
representation of classical statistical models with local 
interactions. It is well known that the one-dimensional 
spin- 1/2 Ising model can be rigorously solved by express- 
ing its partition function as a product of transfer matrix. 
This is the simplest example of tensor network represen- 
tation of classical statistical model. A matrix is a order 
two tensor. A tensor network representation of classi- 
cal statistic model is a higher dimensional extension of 
the one-dimensional matrix product. All statistical mod- 
els with short-range interactions, such as the Ising model 
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and the Pott's model, can be expressed as tensor network 
models. 

There are two approaches to represent the partition 
function of a classical statistical model in a tensor- 
network form. One is to take a duality transformation of 
the model and define the tensor in the dual space. 3 This 
can transform a classical statistical model to a tensor- 
network model in its corresponding dual lattice. The 
order of the local tensor is the coordinate number of 
the dual lattice. This kind of transformation is partic- 
ularly useful if the coordinate number of the dual lat- 
tice is smaller than the coordinate number of the origi- 
nal lattice. For example, a classical model defined on a 
triangular lattice whose coordinate number is 6 can be 
represented as a tensor-network model on a honeycomb 
lattice whose coordinate number is 3. The dual lattice of 
a dice lattice is a Kagome lattice. The square lattice is 
self-dual. 

If the original statistical model includes not only the 
nearest neighboring two-site interactions, but also many- 
site interactions within each constitutional unit cell, for 
example three- or four-site interaction within each pla- 
quette in a square lattice, a tensor-network model can 
still be defined in the dual lattice without enlarging the 
bond dimension of tensor. This is another advantage for 
representing the partition function in the dual space. 

The second approach is to define the tensor in the origi- 
nal lattice by taking singular value decompositions for the 
local bond partitions. This approach can be also imple- 
mented in all kinds of lattices. The order of the tensor is 
equal to the coordinate number of the lattice if there are 
only interactions between neighboring sites. This kind of 
representation is more appropriate for a lattice, such as 
a honeycomb lattice, whose coordinate number is small. 

A. Tensor network representation in the dual 
lattice 

Let us take the 5=1/2 Ising model defined on a trian- 
gular lattice as an example to illustrate how to express its 
partition function as a tensor product in its dual lattice. 
The Hamiltonian of the Ising model is defined by 

H = -jJ2s i S j , (3) 

where Sj takes two values ±1. This model is rigorously 
soluble |39|. 

The partition function of the Ising model is given by 

Z = Tre-P" 

= Tr Y[ eP J( - SlSl+s > Sk+SkSl)/2 ' (4) 
A ijfc 

where Tr is to sum over all spin configurations and the 
product is taken over all small triangles. The factor 1/2 
is introduced in the above exponent because each bond is 
shared by two triangles. This partition function can be 



expressed as a tensor network in terms of bond variables 
through a duality transformation. The dual lattice of a 
triangular lattice is a honeycomb lattice. 

Given a bond on the triangular lattice, let us define 
the bond spin by the the product of two end spins as 

o%j = SiSj. (5) 

This bond spin takes two values a = 1 or —1, correspond- 
ing to a state of parallel or antiparallcl Ising spins. The 
partition function can then be expressed as 

Z = Tr JJ<S((7y - S l S J ) TJ e fiJ ^ + ^ k+a ^' 2 , (6) 

(ij) &ijk 

where Tr is to sum over all S and a spins. 

Since the number of the bonds (Nb on d) is equal to the 
sum of the number of sites (iV s , te ) and the number of tri- 
angles {Ntriangie) , the Nb on d delta functions in the above 
equation can then be equivalcntly written as a product 
of N S ite delta functions defined on all lattice sites and 
Ntriangie delta functions defined on all triangles. On each 
triangle, it is simple to show that the product of the three 
bonds spins is equal to one 

<?ij&jk.Vki = SiSjSjSkSkSi = 1. (7) 

This is the constraint defined on each triangle, indepen- 
dent on the original spin variables. Thus the N S it e -site 
delta functions can be integrated out. The partition func- 
tion then becomes 

Z = Tla II 1 + a V a i kCJki e PJ(viS+<r } k+<7ki)/2_ (g) 

2 

Now if use i to denote the site position of the dual lat- 
tice, then the above partition function can be expressed 
as a standard tensor-network model 

Z = TrJjT Wi) (9) 

i 

where the trace is to sum over all indices and 

T XiViZi = 1 + a *£v* (7 « e m°* i +°y i +°* i )/* (10) 

is a third-order tensor defined on the hexagonal lattice. 
Xi, yi, and Zi are the three integer bond indices of di- 
mension D c = 2 defined on the three bonds emitted from 
site i along the x, y, and z directions, respectively. Each 
bond links two sites. The two bond indices defined from 
the two end points take the same values. 

The above derivation can be extended to other lat- 
tices. A square lattice is self-dual. For the spin-1/2 
Ising model, it is straightforword to show that the tensor- 
network representation of the partition function on the 
square is given by 

Z = T X .y iX l .y'. , (11) 
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FIG. 1: (color online) Schematic representation of the parti- 
tion function of a classical system defined by Eq. ([8| on the 
dual lattice of the triangular lattice, namely the honeycomb 
lattices. 



where 



(12) 

is a fourth-order tensor defined on the dual square lattice. 
Xi, jji , x'i and y\ are the four bonds connecting site i. 

For the honeycomb lattice, its dual lattice is a triangu- 
lar lattice. In the dual space, each site has six neighbors. 
The order of the local tensor is six. Thus there is no 
advantage to carry out the calculation in the dual space. 

The above duality transformation can be extended to 
other statistical models. For example, the q-statc Pott's 
model is defined by the Hamiltonian on a triangular lat- 
tice 



(13) 



(id) 



where Oi = 0, 1, • • • q—1 and S(a—b) is the Kronecker delta 
function, it is straightforward to show that its partition 
function can be also written as a tensor product in the 
dual lattice. But the local tensor is now defined by 



Tr 



8{m.od[x + y + z, q ]) e -PJlHx)+s(v)+Hz)]/^ ( 14 ) 



where the bond variable x (similarly y or z) takes integer 
values from to q — 1. When q = 2, the Potts model is 
equivalent to the spin-1/2 Ising model. 

However, it should be pointed out that Eq. §5§ is not 
a one-to-one mapping. The above derivation can be ap- 
plied to the Ising model without magnetic field. If a 
Zeeman term or a many-site interaction term is added to 
the Hamiltonian, then the transformation Eq. (O can no 
longer be used. Nevertheless, the tensor-network model 
can still be defined in the dual lattice. But the dimension 
of the tensor has to be extended from D c = d (d = 2 for 
the spin-1/2 Ising model) to D c = d 2 = 4 to distinguish 
the 4 spin configurations of Si and Sj. 

To understand this more concretely, let us consider an 
extended spin-1/2 Ising model in an applied magnetic 
field defined on a square lattice 



H = -J^2s i S j -h^2Si + J n ]T 

(ij) i ijkien 



SiSjSkSi, (15) 



Si 



"Oil 



Oik 



FIG. 2: (color online) A unit cell of a square lattice and its 
corresponding dual variables. 



where the last term is to sum over all four-spin inter- 
action terms defined on each plaquette. ijkl E □ means 
that ijkl arc the four-vertex indices on each plaquette, as 
shown in Fig. [2j The partition function of this model can 
still be expressed as a tensor product defined by Eq. (|lljl . 
But the local tensor T x . y . x ' y { is now defined by 



(16) 



where 



2 i^iSi + SiSj + SjSk + SiSk) 
-~ (Si + Sj + Si + S k ) + JaSiSjSkSi. (17) 



If we assume Si = ±1 for the spin-1/2 Ising model, then 
Iij is defined by 

Ui = {Si + 1) + (Sj + l)/2. 

It takes four integer numbers from to 3. The bond 
dimension of the tensor is D c = 4. 



B. Tensor network representation in the original 
lattice 



A classical statistical model with local interactions 
in any spatial dimension can be always expressed as a 
tensor-network model in its original lattice. To do this, 
let us consider a model defined by the following Hamil- 
tonian 



H = J2K(S i ,S j ), 

(ij) 



(18) 



where ( ) denotes summation over nearest neighbors only. 
K(Si, Sj) is the interaction between two local basis states 
Si and Sj . The partition function of this model is given 
by 



where 



Z=T,I[W(Si,Sj), 

is,} m 



W(S i ,S j )=exp[-0K(S i ,S j )] 



(19) 



(20) 



defines a matrix whose row and column indices are Si and 
Sj , respectively. W is not necessary to be symmetric and 
positive denned. 

Given W(Si, Sj), one can take a singular value decom- 
position to decouple it into the following form 

WiSuSj) =J2u(S i J)\iV(S j ,l), (21) 
i 

where U and V are unitary matrices, A is a semi-positive 
diagonal matrix. If we define 

Q a (S,l) = U(S,l)\] /2 , (22) 
Q b (S,l) = V(SJ)\y 2 , (23) 

then W can be reexpressed as 

W = Q a Q b . (24) 

Now let us group all Q's that connect to site i. A local 
tensor can then be defined by tracing out Si from the 
product of these Q's 

T x y z = ^ ' Qctx (Si, x)Qa 2 (Si, y)Qa 3 (Si, z) ■ ■ ■ . (25) 

Si 

The order of the tensor is equal to the number of bonds 
interacting with site i. The bond dimension is equal to 
the site dimension d. This gives a tensor-network repre- 
sentation of the partition function 

y- i I'll (26) 

i 

The two bond indices defined from the two end points 
take the same values. The trace is to sum over all bond 
indices. 

On a bipartite lattice, for example a honeycomb lattice, 
the partition function can be simply expressed as 

Z 11 II TZ.voAv,,.,* (27) 

where the superscripts a and b stand for the two sublat- 
tices of the honeycomb lattice, and 

Ky,* = 5lQa{S t x)Q a (S,v)Q a (S t z), (28) 
s 

T b x , v , z = Y,Qb(S,x)Q b (S,y)Q b (S,z). (29) 

s 

Moreover, if W is positive and symmetric, then V = U 
and Q a = Q b . In this case T = T a = T b and the partition 
funciton can be simply expressed as 

Z Tr J {/,.„.. . (30) 
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(a) (b) (c) 

FIG. 3: Rewiring of a honeycomb lattice by the singular value 
decomposition in the TRG iteration. 



III. SECOND RENORMALIZATION OF 
TENSOR NETWORK MODEL 

A. Overview of TRG 

We first review briefly the TRG methodQ. Let us take 
the model defined by Eq. f|2T[) on a honeycomb lattice as 
an example to show how the method works. 

The TRG is an iterative method. There are two steps 
at each iteration. The first step is to rewire a honey- 
comb lattice to a triangle-honeycomb lattice as shown in 
Fig. [3](b). This is done by transforming a pair of neigh- 
boring tensors, T° and T b , into two new tensors S a and 
S b defined in the rewired lattice, by performing a singu- 
lar value decomposition as schematically shown in Fig. [4] 
The common bond index of T a and T b is first contracted 
to form a matrix M defined by 

m 

The dimension of M is the product of the corresponding 
bond dimensions. The initial dimension of M is D\, The 
singular value decomposition is then applied to decouple 
this matrix into the following form 

Mujk = ^ Uli tn A n Vjk,n, (32) 
n=l 

where U and V are two unitary matrices. A« is a semi- 
positive diagonal matrix arranged in descending order. 
It measures the entanglement between U n and V n . The 
dimension of A is equal to that of M, higher than the 
original bond dimension. To carry out the calculation 
iteratively, one has to truncate the basis space to a man- 
ageable level and retain D cut largest singular values and 
the corresponding vectors. After that, one can define two 
new tensors 

Snii = Uu >n \/ A n , (33) 
S b njk = V jk , n JK~ n . (34) 

The second step is to contract each small triangle in 
the rewiring lattice to define two new local tensors in the 
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FIG. 4: Schematic representation of Eqs. (|31|l and (|32|l . 
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FIG. 5: Schematic representation of Eq. (|35(l . A local tensor 
T a is defined by contracting the three internal bonds on each 
triangle. 

squeezed honeycomb lattice (Fig. O 

rpa \ oct na qa f^l^\ 

1 xyz / j ^xji^ykj J ziki V* 30 7 

where a = a or b. After that, the lattice size is reduced 
by a fact of 3. 

By repeating the above steps iteratively, one can fi- 
nally reach a hexagonal lattice with only six tensors. The 
partition function can then be evaluated by tracing out 
all bond indices of these six tensors, assuming a central 
symmetric boundary (Fig. [S]). 

B. SRG 

The above TRG iteration minimizes the truncation er- 
ror of local matrix M. However, as pointed out in Ref. @ , 
it does not consider the influence of the environment 
which contains all the lattice points excluding the two 
on which M is defined (Fig. [7]) . In real calculation, it is 
the truncation error of the partition function rather than 




FIG. 6: The last six sites in the TRG iteration. A central 
symmetric boundary condition is assumed. 



FIG. 7: Configuration of a system (a) and its corresponding 
environment lattice (b). 

that of the local matrix M that should be minimized. 
This means that the TRG is just a local optimization 
method. 

To optimize the partition function globally, one needs 
to consider the renormalization effect of the environment 
to M. In Ref. @, this is called the second renormaliza- 
tion method of tensor network model (SRG). This SRG 
method, as demonstrated in Ref. Q, can incorporate ef- 
ficiently the renormalization effect of environment and 
improve significantly the TRG method. 

The difference between the TRG and SRG is similar 
to the difference between the conventional Wilson block 
renormalization group method (38j and the DMRGfTTI]. 
In the conventional block renormalization group method, 
it is the block Hamiltonian that is optimized without con- 
sidering the interaction between different blocks. How- 
ever, in the DMRG, the basis states of the system block 
are optimized by fully considering the interplay between 
the system and environment blocks via the reduced den- 
sity matrix. The singular values in the DMRG are the 
coefficients of the Schmidt decomposition of the density 
matrix. They measure the entanglement between system 
and environment blocks. 

This can be understood more clearly by expressing the 
partition function as a product of M and its correspond- 
ing environment matrix M e 

Z = TrMM e , (36) 

where M e is defined by contracting over all bond indices 
in the environment lattice. From this formula, it is clear 
that to reduce the error in Z, one needs to minimize the 
truncation error of MM e , not just that of M. 

We will discuss how to determine the values of M e 
later. Once M e is known, its renormalization to M can 
be done in the following steps: 

1. To take a singular value decomposition for M e 

M!k, U = J2 U knKVlln, (37) 
n 

where U e and V e are two unitary matrices and A e is 
a semi-positive diagonal matrix. This step is taken to 



ensure that the renormalization effect of M e to M can 
be more symmetrically treated and the truncation error 
of the partition function is minimized. 

2. From the above decomposition, one can define a 
new matrix: 

Mn u n 2 = £ (A e ni ) 1/2 Vt tni M Utjk Uf ktna {K 2 ) 1/2 , 
lijk 

(38) 

and represent the partition function as 

Z = TrM. (39) 

This equation means that the error of the partition is 
minimized if the truncation error of M is minimized. 
Now we perform a singular value decomposition for M 

M ni ,n 2 =^Un 1 ,nA„V„ 2 , n , (40) 
n 

again. U and V are two unitary matrices, and A is a 
semi-positive diagonal matrix. According to the least 
square principle, the truncation error of M is minimized 
if the D cut largest singular values of A are retained in the 
truncation. 

3. By substituting the truncated M back to Eq. [33) 
one can represent M as 

Mu t jk = y^Vij, ni ( A nJ 1/2 M ni>n2 (A^ J 1/2 Uj k 7l2 . 

(41) 

It can be further expressed as a product of two tensors 
M H , jk « J2 S Z,Ak, ( 42 ) 

n=l 

where 

Sin = E V «.nx (Kj'^U^.n (A n ) ^ , (43) 

SL'fc = £^,n 2 ( A «J" 1/2 K 2) „ (A„) 1/2 , (44) 

712 

are the two tensors defined in the rewired lattice. 

Then one can follow the steps of TRG to update ten- 
sors T a and T b in the squeezed lattice by taking the 
coarse grain decimation of S a and S b . This completes a 
full cycle of SRG iteration. By repeating this procedure, 
one can finally obtain the value of partition function in 
the thermodynamic limit. 

C. Determination of the environment tensor M e 

In the DMRG calculation, a superblock can be sepa- 
rated into a system block and an environment block by 
cutting one bond. Thus the environment can be read- 
ily identified and integrated out. However, in the TRG 
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A=A 

FIG. 8: Schematic representation of the iterative renormal- 
ization to the singular values A of M. M is obtained from 
Eq. p8[) by taking a mean-field approximated M e defined by 
Eq. (3S)| . 

method, it is highly nontrivial to handle the environment. 
To separate the environment from the system block, one 
needs to cut four bonds. The system contains only two 
sites and the environment contains all rest of sites. In 
this case, it is even more difficult to evaluate the con- 
tribution of the environment than the partition function 
itself in the TRG. 

The key step in the SRG is to calculate the environ- 
ment tensor M e . In this respect, two approaches can 
be used. One is to take a mean-field approximation (or 
cavity approximation) to account for the environment 
contribution. [40j This is a cheap but less accurate ap- 
proach. It is based on an intuitive interpretation to the 
singular values A of M. The other is a more accurate 
approach. It is to evaluate M e directly from the environ- 
ment lattice. 

Let us first consider the "mean field" approach [40j]. As 
mentioned before, the singular bond vector A„ of M de- 
fined by Eq. (|32|) is a measure of the entanglement be- 
tween the corresponding basis states U n and V n . It can 
be also regarded as a measure of the interaction between 
the two-end basis tensors U and V linked by this bond. 
If we assume that this singular vector also measures the 
entanglement between the system and the environment 
on all four dangling bonds (i, j, k, I) connecting these two 
subsystems, then the environment matrix M e (up to an 
irrelevant prefactor) is approximately given by 

Mf hjk » y/ AiAiAjAk. (45) 

The reason we use A 1 / 2 instead of other exponent of A as 
the weighing factor from the environment is because half 
of A can be associated to the system and the other half to 
the environment. This is a simple but crude approxima- 
tion. However, as will be shown later, it does reveal the 
importance of the renormalization effect of environment 
in the optimization of TRG. 

To apply the above approximation, one needs to first 
calculate the singular vector A of M, and then substi- 
tutes it into Eq. (gSJ) to find M e . From Eqs. (J3SJ) and 
(j4"0)) . one can find a new singular vector A. This A incor- 
porates partial contribution from the environment. To 
treat more accurately the environment contribution, one 
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FIG. 9: Configurations of a finite honeycomb lattice with 6, 
10, 16, or 24 sites. The corresponding environment lattice, by 
excluding the two red vertices, contains 4, 8, 14 and 22 sites, 
respectively. 

can replace A in Eq. (f45|) with A and repeat the above 
calculation itcrativcly. A graphic representation of this 
renormalization procedure is shown in Fig. [5J 

In the above renormalization procedure, the truncation 
error of the partition function can be iteratively reduced. 
The iteration can be terminated when the truncation er- 
ror is less than a desired value. In most of calculations, 
we find that two to three iterations are enough. However, 
at the critical point, more iterations are generally needed. 
The precision of the truncation error is not necessary to 
be set too high. If too many iterations are taken, some 
of the smallest singular values may become smaller than 
the machine error. In this case, further iterations will 
increase rather than reduce the error in the final result. 

One can also treat self-consistcntly the above mean- 
field approximation, namely to require the bond vector A 
used by M e in Eq. (|45|) to be the singular values defined 
in Eq. (|40|) . This self-consistent approximation implies 
that the system is scaling invariant. It might be a good 
approximation at the critical point. 

A more accurate and reliable approach is to evalu- 
ate M e from the environment lattice without taking the 
above mean-field approximation. For small lattice sys- 
tems, for example, with only 6, 10, 16, and 24 sites as 
shown in Fig. [9] (the corresponding numbers of sites in 
the environments are 4, 8, 14, and 22, respectively), M e 
can be rigourously evaluated by contracting over all bond 
indices in the environment. 

Fig. [10] shows how the relative error of the free energy, 

varies with the environment size for the Ising model on 
a triangular lattice. f(T) is the calculated free energy. 




T 



FIG. 10: Relative errors of the free energy for the Ising model 
on a triangular lattice obtained by considering the second 
renormalization effect from four finite environment lattices 
which contains 4, 8, 14, and 22 sites, respectively. The con- 
figurations of these environments are shown in Fig. [9] The 
TRG result is also shown for comparison. 

fex(T) is the exact free energy of the Ising model on a 
triangular lattice derived by Wannier(39l|. The critical 
temperature of a paramagnetic to ferromagnetic phase 
transition is T c = 4/ In 3. As expected, the relative error 
of the free energy decreases with the increase of environ- 
ment lattice. Thus the renormalization effect of M e to M 
can be more and more accurately described by increasing 
the environment size. 

For a larger environment lattice, M e cannot be evalu- 
ated by directly tracing out all bond indices. An accurate 
and efficient approach we proposed in Ref. @ is to use the 
TRG to calculate M e . This can be achieved by perform- 
ing a forward iteration, followed by a backward iteration, 
described below: 

Forward iteration: 

This is similar as in the standard TRG calculation. We 
apply the TRG to coarse grain the environment. At each 
step, one first constructs the M-matrix from Eq. pip by 
tracing out the common bond of two neighboring sites 
on which the M-matrix is defined. Then a singular value 
decomposition for this matrix is taken to find the corre- 
sponding S a ' n and S b - n tensors using Eqs. (fBTj) and ([62]). 
where the superscript n denotes the n'th TRG iteration. 

Fig. [TTJ shows how the environment lattice changes at 
the n'th TRG iteration. For a given environment lattice 
as shown in Fig. fTTT a). a rewired lattice whose configura- 
tion is given by Fig. Illf b) is derived by the singular value 
decomposition. A decimated lattice, shown in Fig.fTlTc'). 
is then obtained by contracting over the three internal 
bonds on each triangle and replacing it by a single lat- 
tice point. This new configuration of environment can 
be separated into two parts as shown in Fig. [TTTd) and 
(e). After rotating by 90 degrees clockwise, Fig. [TTTe) 
looks exactly the same as the original environment lat- 
tice as shown in Fig. ITlT a) , except that the size of this 
reduced environment lattice is by a factor 3 smaller than 
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FIG. 11: One step of the forward iteration of the environment 
lattice. 

its original one. 

If wc use A/(" -1 ) to denote the environment ma- 
trix corresponding to the environment lattice shown in 
Fig- mT a). then from Fig. [TTJ it is straightforward to 
show that it satisfies the following recursion formula 

A/r( n_1 ) _ \^ UW C r»,HnlI,Ilni,Il Q b,U ( . »\ 

m lijk — Z-~j V i r j'k' i'pl j'ip h'qj I'kq 5 I 4 ' I 

I'i'j'k' pq 

where AI e = M^ ' is the final result to be found. 

The above steps can be repeated until M e for a suf- 
ficiently large environment lattice can be accurately de- 
termined. If wc assume that the above iteration is termi- 
nated at the iV'th step so that only 4 sites are left in the 
reduced environment lattice as shown in Fig. [9ja), then 
can be determined by the following formula 

abed 

where the center-symmetric boundary conditions are as- 
sumed. T a (a = a, b) is given by 

rpa \ qOt,N qCt,N qOl,N (AQ\ 

- L xyz / j °xii ykj ° zik ' 
ijk 

Backward iteration: 

In Eq. (|47|) . M^™ -1 ) is determined by the environment 
matrix M^ n > which is to be evaluated from the next iter- 
ation. Thus M^ n ~ x ' cannot be directly determined from 
the above forward iteration. However, we can do a back- 
ward iteration from M^ N ' to find the environment matrix 
AI e = M(°)) using Eq. (|47|) . 

From the above forward-backward iteration, one can 
calculate accurately the environment tensor M e . In prac- 
tical calculations, we find that there is no need to use 
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FIG. 12: (color online) Comparison of the relative error of the 
free energy for the Ising model on triangular lattices obtained 
using TRG (red), the mean-field approximated SRG (blue), 
and the SRG (black) methods with D cu t = 24, respectively. 
The critical temperature is T c = 4/ In 3. 

a large N, since the environment tensor M e converges 
quickly with the increase of environment lattice. 

Fig. [12] compares the relative errors of the free energy 
for the Ising model on a triangular lattice obtained using 
the TRG, the mean- field approximated SRG, and the 
SRG methods, respectively. We find that the mean-field 
approach of SRG can improve dramatically the accuracy 
of the results. This shows that our intuitive assumption 
for the mean-field character of A has indeed caught the 
main feature of the entanglement between the system and 
environment lattices. As expected, the improvement of 
SRG is even more impressive. It improves the accuracy 
for more than five orders of magnitude away from the 
critical point and for more than two orders of magnitude 
at the critical point by keeping only D cut = 24 states. 
The accuracy can be further improved if the environment 
is determined self-consistcntly by SRG which, however, 
will increase the time cost. 

Furthermore, we find that the improvement of the 
SRG over the TRG becomes more and more pronounced 
with increasing D. Fig. [13] shows how the relative er- 
rors change with D cut for the Ising model on triangular 
lattices. 



D. Other lattices 

The above SRG approach can be extended to apply to 
other kinds of two dimensional lattices, such as square or 
Kagome lattices. 

Let us first consider a square lattice. Unlike in a honey- 
comb lattice, the M-matrix is now the T-matrix defined 
at each lattice site. The TRG iterations on square lattices 
can be done using the approach introduced in Ref. [3| • To 
evaluate the environment matrix A4 e at each TRG iter- 
ation, one can still use the forward-backward iteration 
scheme of SRG introduced on a honeycomb lattice. To 
do this, we need to first convert a square lattice to a 
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FIG. 13: (color online) The relative error of the free energy 
as a function of the truncation dimension D cu t for the Ising 
model on triangular lattices obtained using the TRG (black) 
and SRG (blue), respectively. T = 3.2 



honeycomb lattice. This can be done by singular decom- 
posing the fourth order T-tcnsor defined at each square 
lattice point into two separated third order tensors along 
a diagonal direction. A graphical representation of this 
decomposition of the square lattice is shown in Fig. 1141 

After this conversion, the environment has the same 
configuration as that for a honeycomb lattice shown in 
Fig.fTlTa). The only difference is that in the present case 
the dashed bonds in Fig. [T47 b) have higher dimensions 
than the solid bonds. Therefore, the forward-backward 
iteration scheme previously introduced can be used to 
calculate M e for a square lattice. 

Fig. [15] compares the relative errors of free energy for 
the Ising model on the square lattice obtained by the 
TRG with those obtained by the SRG. Similar as for 
the honeycomb lattice, we find that the SRG is much 
more accurate than the TRG. It improves the accuracy 
for more than three orders of magnitude away from the 
critical point and for more than one order of magnitude 
at the critical point over the TRG by keeping D cut = 24 
states. 

The SRG scheme introduced on the honeycomb lat- 
tice can be also applied to a tensor-network model de- 
fined on a Kagome lattice. A Kagome lattice can be 




FIG. 14: To convert a square lattice (a) to a honeycomb lat- 
tice (b) by singular decomposing the fourth order tensor T at 
each point into two separated third order tensors along the 
direction indicated by the red dashed lines. 
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FIG. 15: Comparison of the relative errors of the free energy 
for the Ising model on a square lattice obtained by the SRG 
with those obtained by the TRG. D = 24 

converted to a hexagonal lattice in the following two 
steps: First, we singular decompose each local tensor 
to rewire the Kagome lattice (Fig. [T6a ) to a triangular- 
hexagonal lattice (Fig. ITfjb); Second, we contract all tri- 
angles in Fig. [Tfjr b) by tracing out all internal bonds on 
each triangle. This yields a honeycomb lattice as shown 
in Fig. [ToTc). Then one can calculate all thermodynamic 
quantities just using the method developed for the hon- 
eycomb lattice. 

IV. QUANTUM TENSOR NETWORK STATES 

As mentioned before, the wavef unction generated by 
the DMRG is a matrix product state. It can be shown 
that any wavefunction in one dimension can be faithfully 
expressed as a matrix product. In two or higher dimen- 
sions, a matrix product wavefunction is no longer a good 
description of the ground state, since the minimal ma- 
trix dimension that is needed for accurately describing 
the state grows exponentially with the lattice size. 

A tensor network state is a natural extension of the 
one-dimensional matrix product state in two or higher 
dimensions. It captures the main feature of entanglement 




FIG. 16: To convert a Kagome lattice (a) to a honeycomb 
lattice (c) by rewiring with singular value decompositions and 
contraction. 
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between different parts of wavcf unction. For a quantum 
spin system, for example, a tensor network wavefunction 
can be expressed as: 



(50) 



where T%. mZi ...[mi] is the tensor defined at site i. The 
order of this tensor is the number of spins interacting 
with Si. If there are only nearest- neighbor interactions, 
the rank of T l [mj\ is just the coordinate number of the 
lattice, which is 3 for a honeycomb lattice, 4 for a square 
or Kagome lattice, and 6 for a triangular lattice. The 
subscripts (xi, yi, Zi ■ ■ •) are bond indices, to.; is a set 
of local basis states. The trace is to sum over all basis 
configurations and over all bond indices. The local tensor 
T 1 is generally site dependent. However, if the system is 
translational invariant, then T l can be site independent. 

To determine the local tensor T 4 [m,], a commonly 
adopted approach is to take all the tensor elements as 
variational parameters and to determine them by varia- 
tionally minimizing the ground state energy 



E 



(51) 



However, the number of variational parameters that can 
be efficiently handled by the minimization is only a few 
hundreds. This limits the bond dimension that can be 
handled to be generally less than 5 on a square lat- 
tice. Moreover, the accuracy of the wavefunction such 
obtained is very low. 

In Ref. [H, we proposed an iterative projection method 
to determine the tensor-network wavefunction. This is 
a more efficient and accurate method. Below we take 
the S=l/2 Heisenberg model on a honeycomb lattice to 
demonstrate how the method works. The model Hamil- 
tonian is defined by 



H 



Hi 



(ij) 

J Si ' Sj , 



(52) 
(53) 



where (ij) stands for summation over nearest neighboring 
sites. 

The honeycomb lattice is a bipartite lattice. It can be 
divided into two sublattices, denoted as a and 6, respec- 
tively. A translation invariant tensor-network wavefunc- 
tion can in general be represented as: 

\ip)=Ti X Xi X yi X Zt A XiViZi [Tn>i\B Xj y jZj [uij] \miTUj) 

i£a : j£b 

(54) 

m,j\ are the tensors defined 
on the two sublattices, respectively. A a (a = x,y,z) is 
a positive bond vector of dimension D, defined on the 
bond emitted from site i along the a direction, to, is a 
local spin basis set of dimension d. The trace is to sum 
over all spin configurations and over all bond indices. 



where A XlVlZt [to, ; ] and H x v • [ 



In Eq. (|54p . we introduce explicitly the bond vectors 
A in the tensor network wavefunction. These A, as will 
be discussed below, measure approximately the entan- 
glement between two neighboring sites. Since each bond 
connects with two tensors, one can regard each A as a 
product of two A 1 / 2 and associate each of them to one of 
the tensors. In other words, one can take A 1 / 2 as a mean 
field measure of entanglement acting from one tensor to 
another. In the calculation of the tensor network state by 
projection, the renormalization effect from the environ- 
ment tensors need to be considered. This renormaliza- 
tion effect, similar as in the application of the SRG, can 
be taken into account approximately by these effective 
entanglement fields. This approximate treatment of the 
renormalization effect from the environment, as will be 
discussed below, works very well. It provides an accurate 
and efficient method for determining the tensor-network 
wavefunction. 

In the projection method, the ground state wavefunc- 
tion is determined by applying the projection operator 
exp(— tH) to an arbitrary initial state which is not 
orthogonal to the true ground state. In the limit r — > oo, 
the resulting wavefunction will converge to the ground 
state of H. However, this projection cannot be done in 
a single step since the terms in H defined by Eq. (f5"3")) 
do not commute with each other. Instead, one needs to 
use a small r and apply this projection operator to \^J) 
itcratively for sufficiently many times. 

Let us start by dividing the Hamiltonian into three 
parts 

H = H x + H y + H z , 
H a = J~] Hj ti+a (a = x,y,z). 

H a (a = x, y, z) contains all the interaction terms along 
the a-direction only. These terms commute with each 
other. For small r, one can use the Trotter-Suzuki for- 
mula to decouple approximately exp(— tH) into a prod- 
uct of three terms 



-tH 



+ o(r 2 ). 



(55) 



A projection of H can then be readily performed using 
cxp(— TH a ) (a = x,y,z) in three steps. 

In the first step, the projection is done with H x . As H x 
contains only the interaction terms between two neigh- 
boring spins connected by horizontal bonds, this projec- 
tion generates a new wavefunction 



Tr 



I*) 

n e 

i^a,j— i+x mifrhj 



m'im'jle HijT \mimj 



K i \Az t A Xz y iZi [m i \B XjyjZj [m j }\m' l rn' j ) 1 (56) 

where any two tensors linked by a horizontal bond are 
mixed by the matrix elements of the local projection op- 
erator. In order to perform the projection for the next 
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step, one needs to separate these two tensors so that the 
wavefunction can return back to its original form. 

To do this, let us define a {D 2 d) x (D 2 d) matrix us- 
ing the two tensors on a horizontal bond and the bond 
vectors connected to these tensors 

<? / / 

^yiZiTn,^ ,yj Zj rrij 
rriimj x 

[mi]X x B X y jZj [mj]X yj X zr (57) 

In this definition, the four bonds that connect to the 
environment, i.e. yi, Zi, yj and Zj, are weighted by the 
corresponding bond vector A, rather than A 1 / 2 . These 
extra A 1 / 2 ^ arc included on these four bonds to mimic 
the renormalization effect from the environment. This, 
similar as in our Poor Man's treatment of SRG, is a mean- 
field- type treatment to the environment. It is a key step 
in the application of this projection method. Without 
considering this renormalization effect, the iteration does 
not even converge. 

We then take the singular value decomposition to de- 
compose the above 5* matrix into a product of two tensors 
at sites i and j, respectively 

Sy i Zinii^yjZjmj — ^ ^ UyjZjrnj ,xXxVyj Zjnij ,x (^^) 

X 

where U and V are two unitary matrices and \ x is a 
positive diagonal matrix of dimension D 2 d. X x is the 
singular value matrix of S. It measures the entanglement 
between U and V tensors. 

To process the iteration, one has to truncate the basis 
space by keeping only the D largest singular values of 
X x . The truncated X x is then set as the new X x for the 
next iteration. After this, tensors A and B are updated 
as follows 

A xyiZi [mi\ = Xy^X^Uy iZzrnuX , (59) 
B xyjZj [mj} = X~jX^V VjZjmjtX . (60) 

This completes the projection for all horizontal bonds. 
The projection for the bonds along the y- or z-direction 
can be done in the same way. By repeating this iteration 
procedure, an accurate ground state wave function can be 
finally determined. In practical calculation, we find that 
the projection converges more efficiently if the second- 
order Trotter-Suzuki decomposition formula is used to 
decouple the projection operators along different direc- 
tions. 

In Eq. (|57| . the renormalization effect of the environ- 
ment is taken into account by a mean-field approxima- 
tion. This approximate treatment does not minimize the 
truncation error at each step of projection. However, this 
does not affect the converging speed of the wavefunction, 
since the Trotter and truncation errors do not accumulate 
in the iteration of projection. 




Iteration Number 

FIG. 17: (color online) Upper panel: Comparison of the 
ground state energy of the ID Heisenberg model as a func- 
tion of iteration number obtained using the projection method 
with and without taking canonical transformation for the 
matrix-product wavefunctions, E can and Emf- Lower panel: 
the energy difference Emf — Ecan- The bond dimension 
D — 10. 



To see how well this projection method works, let us 
consider a one-dimensional system. In one dimension, 
one can do a canonical transformation to convert a ma- 
trix product wavefunction into a canonical form, in which 
the left and right subblocks are strictly orthogonal and 
A 2 are the eigenvalues of the reduced density matrix |37}. 
If this canonical transformation is applied to the matrix 
product state before each projection, then the truncation 
error is minimized and the wavefunction such obtained is 
the most accurate one. 

In the upper panel of Fig. [T71 we have compared the 
converging speed for the ground state energy obtained 
by taking the canonical transformation, E can , with that 
obtained by simply using the mean-field approximation, 
Emf, starting from a randomly generated wavefunction 
for the one-dimensional Heisenberg model with the bond 
dimension D = 10. As can be seen from the figure, both 
E can and Emf converge to the true ground state energy 
of D cut = 10 quickly. It is impossible to see the difference 
between E can and Emf on the energy scale of the figure. 
In the lower panel of the figure, we show the energy dif- 
ference Emf — E can as a function of the iteration num- 
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ber. At the beginning of the iteration, Emf is slightly 
higher than E can , indicating that the wavefunction ob- 
tained by the canonical approach is indeed more accurate, 
as expected. However, the difference between these two 
energies approaches to zero with increasing number of it- 
erations, indicating that the wavefunction obtained using 
the mean- field approach does converge to the true ground 
state wavefunction. The application of the mean-field 
method does not need to perform the canonical transfor- 
mation at each step of projection. Thus in many cases it 
can be more efficient than the canonical approach in the 
evaluation of the wavefunction. 

However, it should be pointed out that the above 
mean-field treatment to the environment cannot be sim- 
ply applied to the calculation of time-dependent or ther- 
modynamic quantities with the above projection method. 
This is because in the calculation of time-dependent or 
thermodynamic quantities, both the Trotter and trun- 
cation errors will accumulate with the step of iterations 
and the long-time (or imaginary time in the calculation of 
thermodynamic quantities) results may not converge. In 
this case, a more rigorous treatment to the renormaliza- 
tion effect from the environment should be considered to 
minimize the accumulated Trotter and truncation errors. 

In practical calculation, the symmetry of the Hamil- 
tonian can be used to block diagonalize tensors A and 
B in Eq. (|54"|) . This can improve greatly the efficiency 
of the calculation and allow larger bond dimension to be 
handled. For example, the conservation of the spin com- 
ponent along the z-axis, S z , can be readily enforced by 
requiring the bond spins to satisfying the equation 

S[xi] + S[yi] + S[zi] = m l; (61) 

for the tensor on the a-sublatticc, and the equation 

-S[x j ]-S[y j ]-S[z j ]=m j . (62) 

for the tensor on the 6-sublattice. In the above equations, 
S[a] is the spin quantum number of the bond index a. 

Graphically, the above spin conservation law can be 
represented by an arrow rule as shown in Fig. 1181 We 
put an arrow on each bond to represent the sign factor 
for the virtual bond spin. An inwards/outwards arrow 
takes a positive/negative sign. The conservation law is 
defined at each vertex, given by the equality that the 
sum of the bond spin times its arrow sign is equal to 
the physical spin at the vertex. This arrow rule can be 
extended to apply to all kinds of lattices, no matter they 
are bipartite or non-bipartite. 

At each projection of a pair of tensors, for example 
the two tensors linked by a horizontal bond, A XilVilZi and 
Bx^yjZj {xj = Xi), the conservation of total S z is satisfied 
if ' 3 3 



S[yi] + S[zi] - m, = S[yj] + S[z 3 } + m } 



(63) 

From Eqs. (|59|) and ([60|) . it is simple to show that the 
sum of S z on the whole lattice is zero: 





FIG. 18: (color online) Graphical representation of the spin 
conservation law denned by Eqs. Q6ip an d (|62p . The bond 
spin takes a positive/negative sign if the bond arrow points 
towards/outwords the vertex. 



Thus the tensor-network wavefunction satisfying the 
above conservation conditions has total S z = 0. One 
can also construct a tensor-network wavefunction with 
nonzero total S z by modifying Eq. (|59|) or (|60|) on some 
of the lattice points. But in this case, the wavefunction 
is no longer translation invariant. 

It is straightforward to extend the above projection 
method to other quantum lattice models with short range 
interactions either in a honeycomb or other kinds of lat- 
tices in two or higher dimensions. At each step of projec- 
tion, a singular value decomposition of a D z ~ 1 d x D z ~ 1 d 
matrix needs to be calculated, where z is the coordinate 
number of the lattice. The largest bond dimension D 
that can be studied depends on the coordinate number 
z. The larger z, the smaller D that can be handled. 



V. CALCULATION OF PHYSICAL 
OBSERVABLES 

Now let us consider how to calculate the expectation 
values of a physical operator O 



(0) = 



(*|Q|tt) 



(65) 



For a given tensor-network wavefunction |\&), it is simple 
to show that both the denominator and numerator can be 
expressed in the form of tensor products. In particular, 
the denominator has a similar form as for the partition 
function of a classical statistical model: 

Tr JJ V.v yz's yz- (66) 

where T a and T b are the local tensors on the two sublat- 
tices defined by 



(67) 



0. 



(64) 



T x t Y iZi = E^wH^H, (68) 



Xi = (x 4 ,x-), Yi = {yi.y'i), and = (z^z'A. The bond 
dimension of these local tensors is D r — D 2 . 
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FIG. 19: (color online) The SRG result of the ground state 
energy as a function of the truncation dimension D cut for the 
Heisenberg model on a honeycomb lattice. D is the bond 
dimension of the wavefunction. 



The numerator differs from the denominator only in 
the definition of the local tensors on the sites where the 
physical operators are defined. For example, in the eval- 
uation of the ground state energy, one needs to calculate 
the expectation value of H\2'. 



(69) 



This can be calculated from the sum of the above three 
terms on the right hand side. The corresponding local 
tensors on sites 1 and 2 are defined by 

T x-r i z i = ^ A xiViZj M A^g/ [m'\ (m'\ Si lCt |m) , 



T 



XiYiZi 



^ BxiViZi M B x' i y' i z' f W] (m'\S 2 , a \m) 



The definition of local tensors on the other sites is un- 
changed. For the Heisenberg model on the honeycomb 
lattice, the ground state energy per-site is given by 
3<#i2>/2. 

The spin conservation law defined by Eqs. (|61|) and 
(p32"j) can be also implemented to tensors T a and T b . If 
we define 



S[Xi] = S[x t ] - Slx'j] 



(70) 



and similarly for S[Yi] and S 1 ^], then the spin conser- 
vation for T a or T b is then given by 



S[Xi] + S[YJ + S[Zi] = 0. 



(71) 



The expectation value (|65|) can be calculated using the 
SRG method introduced in Sec. HTH Fig. Ql shows 
how the ground state energy varies with the truncation 
dimension D cut for the Heisenberg model on a honey- 
comb lattice. For the three cases shown in the figure, 
the ground state energy shows a small variance at the 



FIG. 20: (color online) The ground state energy of the Heisen- 
berg model on a honeycomb lattice as a function of the bond 
dimension D obtained by the SRG with D cut = 130. 



fifth decimal when D cut is above 80. This variance re- 
sults from the truncation error in the SRG calculation. 
The ground state energy does not decrease monotonically 
with D cut , because SRG is not a variational approach. 

Fig. [20] shows the ground state energy of the Heisen- 
berg model as a function of the bond dimension D ob- 
tained using the SRG. The result converges quickly with 
increasing D. This is because the ground state energy is 
determined by the local correlation function. It is insen- 
sitive to the long wavelength fluctuation of the wavefunc- 
tion. The converged ground state energy of D = 16 is 
-0.54440 by keeping D cut = 130 states. It agrees with the 
most recent Monte Carlo result E = -0.54455(20). [Hj It 
is also consistent with the spin wave[42j (-0.5489) as well 
as the series expansion [43j (-0.5443) results. 

For the Hamiltonian, or any other physical operators 
which commute with H, the ground state |^) can be a 
common eigenfunction of these variables. The expecta- 
tion values of these conserving physical variables can be 
also evaluated by the following equation 



(0) = 



($1*) 



(72) 



where | $) is an arbitrary wavefunction that is not orthog- 
onal to | , F). An advantage for evaluating the expectation 
value using this formula is that the bond dimension of |4>) 
can be much smaller than |\&). This can reduce the bond 
dimension D c for the tensors used in Eq. (|66[) and allow 
a tensor-network wavefunction with a relatively larger 
bond dimension to be studied. 

For the Heisenberg model on a honeycomb lattice, the 
ground state is spontaneous symmetry broken. It pos- 
sesses a long range antiferromagnetic order with a finite 
staggered magnetization. The staggered magnetization 
measures the long range spin-spin correlation functions. 
In an applied staggered magnetic field, the spin wave ex- 
citation is gapped. The spin-spin correlation function, 
excluding a constant long range term, decays exponen- 
tially with their distance. However, in the absence of an 



16 




s 



FIG. 21: (color online) The staggered magnetization as a 
function of D for the Heisenberg model on a honeycomb lat- 
tice. 



FIG. 22: (color online) The field dependence of the staggered 
magnetization M s for the Heisenberg model on a honeycomb 

1/2 

lattice. The inset shows the h 3 dependence of M s . 



applied staggered magnetic field, the spin wave excita- 
tion is a gapless Goldstone mode. The low energy, or 
long wavelength, spin fluctuation is strong in a Honey- 
comb lattice. These low energy fluctuations can affect 
strongly the behavior of the staggered magnetization in 
the low field limit. 

In the ground state of which the SU(2) spin rotation 
symmetry is broken by applying a staggered magnetic 
field h s (assuming along the z-axis), the staggered mag- 
netization can be evaluated from the formula 



1 



(<Si,z — S\ 



2,z) 



(73) 



Fig. [22] shows the SRG results [44J for the staggered 
magnetization M s as a function of the bond dimension 
D in the zero field limit h s — > 0. In contrast to the 
ground state energy, M s converges slowly with increas- 
ing D. When D = 16, we find that M s w 0.3098, which is 



higher than the recent quantum Monte Carlo result [41 
M s « 0.2681(8). It is also higher than the result obtained 
from the spin wave theory [42j|. M s = 0.24, or the series 
expansion [4J], M s — 0.27. This difference is due to the 
quantum fluctuation in the ground state. In the ground 
state, the spin excitation is gapless and the spin-spin cor- 
relation is long ranged. However, for a tensor-network 
wavefunction, this long-range spin-spin correlation is ter- 
minated by the finite bond dimension. This is a drawback 
of the tensor product wavefunction in the studying of a 
gapless state. This kind of error can be reduced by in- 
creasing the bond dimension of the tensor product state 
or by adopting other kinds of tensor-network wavefunc- 
tions. 

To understand this more clearly, we compare the SRG 
result for the staggered magnetization as a function of 
an applied staggered field h s with that obtained from 
the spin-wave theory in Fig. [22] At high field, the two 
results agree qualitatively with each other. This is be- 
cause in a finite field, the spin wave excitation is gapped 
and the entanglement entropy satisfies the area law. In 



this case, the tensor-network wavefunction is a good ap- 
proximation to the true ground state. However, in the 
limit h s — > 0, the spin wave excitation becomes gapless 
and the critical spin fluctuation becomes important. In 
particular, as shown in the inset of Fig. [22] the staggered 
magnetization obtained by the spin wave theory varies 
as \fh~ s in the low field limit. This \fh~ s dependence of 
the staggered magnetization is due to the low-energy (or 
long wavelength) spin excitations. The long wavelength 
correlation is not included in the tensor-network approx- 
imation of the ground state. It leads to the error in the 
SRG result of the zero-field staggered magnetization for 
finite D. 

The staggered magnetization in the zero field limit can 
nevertheless be more accurately estimated by extrapolat- 
ing the SRG result to the limit D — > 00. Fig. [23] shows 
the fourth order polynomial fit to the SRG results of 
staggered magnetization as a function of 1/D. The ex- 
trapolated staggered magnetization in the limit 1 /D — > 
is about 0.285, in agreement qualitatively with the quan- 
tum Monte Carlo result, M s w 0.2681(8), as well as the 
series expansion one, M s = 0.27. 



VI. SUMMARY 

In this paper, we have discussed the tensor-network 
representation of classical statistical models, and given 
a comprehensive introduction to the iterative projection 
and the SRG methods. A classical statistical model with 
local interactions can be represented as a tensor-network 
model either in its original lattice or in its dual lattice. 
The order of local tensor is equal to the coordinate num- 
ber of the lattice on which the tensor-network model is 
defined. In practical calculation, one should choose a 
tensor-network representation in which the order of local 
tensors is the smallest. For example, in the study of a 
classical model in the honeycomb/triangular lattice, the 
tensor- network model defined in the original/dual lattice 
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FIG. 23: (color online) The staggered magnetization as a 
function of 1/D for the Heisenberg model on a honeycomb 
lattice. 

should be used. 

The projection method introduced in Sec. lIVl is an effi- 
cient and accurate tool for evaluating the tensor-network 
wavefunction for a quantum lattice model. It allows a 
tensor-network state with a large bond dimension, for ex- 
ample D = 70 in a honeycomb lattice, to be determined. 
In the projection, the renormalization effect of the envi- 
ronment is taken into account by a mean-field treatment 
of the bond vector. This reduces significantly the trun- 
cation error at each step of projection and enable the 
wavefunction to converge fast with the increase of itera- 
tions. As both the Trotter and truncation errors do not 
accumulate in the iteration of projection, the accuracy of 
the wavefunction can be well controlled by adjusting the 
Trotter parameter r. 

The SRG is an accurate numerical method for evaluat- 
ing thermodynamic properties of classical tensor-network 
models as well as the expectation values of quantum 
tensor-network states. It generalizes the TRG method 
of Levin and NaveQ to account for the renormalization 
effect from the environment in the decomposition of lo- 
cal tensors. This method reduces dramatically the trun- 
cation error and improves significantly the accuracy of 
TRG. The concept of SRG is ubiquitous. The key idea is 
to maximize the entanglement between the system and 
environment in the basis truncation. It can be applied 
to the tensor network models. It can be also extended to 
apply to other physical problems where the system can 
be divided into two parts and the interplay between them 
is important. 

Both the projection method and the SRG can be ap- 
plied to a finite lattice system. They can also be applied 



to a system without any translation or rotation symme- 
try. In this case, one has to treat each tensor indepen- 
dently. The cost (both the CPU time and memory space) 
scales linearly with the lattice size. 

The tensor network wavefunction provides a good de- 
scription for the ground states of two-dimensional quan- 
tum lattice models. For the spin-1/2 Heisenberg model 
on a honeycomb lattice, the ground state energy obtained 
from the tensor network states converges fast with the 
bond dimension D. Our SRG results for D < 16 al- 
ready reach the accuracy of the recent quantum Monte 
Carlo calculation. By further increasing D or reducing 
the truncation error in the SRG calculation, we believe 
that more accurate results for the ground state energy 
will be obtained in near future. 

The tensor network state satisfies the entanglement 
area law. It captures the key feature of short range cor- 
relations. The correlation function between any two lo- 
cal operators in a tensor network state is always short 
ranged. In other words, all low energy excitations are 
gapped in a tensor network state. However, long range 
correlations are not correctly described by a tensor net- 
work state with finite bond dimension. This leads to 
a relatively large error in the determination of physical 
quantities, such as the staggered magnetization in the 
honeycomb Heisenberg model, whose values are governed 
by gapless low-lying excitations. To resolve this problem, 
one has either to increase the bond dimension or to adopt 
a new type of tensor network wavefunction, such as the 
multiscale entangled tensor network state [45j. 

The iterative projection and SRG methods introduced 
in this paper improve significantly the accuracy and ef- 
ficiency in the study of tensor-network states/models. 
These methods can be used to investigate the ground 
state properties of interacting fermions[46j or frustrated 
quantum spin models. It also has the potential to be 
extended to probe thermodynamic as well as dynamic 
properties of quantum lattice models in two dimensions. 
The application of these methods is still in its early 
stage. [43,I!H Further development of these methods may 
enlighten the route for the exploration of new numerical 
renormalization group methods, and lead to the solution 
of a number of problems that are difficult to be solved by 
other methods. 
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